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Abstract 


The  purpose  of  this  work  was  twofold:  (a)  To  develop  numerical  or  analytical  methods 
capable  of  dealing  with  the  combination  of  three-dimensionality  and  strong  anisotropy 
that  occurs  in  the  plume  region  near  the  exit  of  a  Hall  thruster,  and  (b)  To  explore  the  use 
of  fluid-based  models  in  plume  computations,  so  as  to  overcome  the  granularity 
associated  with  normal  particle-based  approaches.  In  area  (a)  we  have  carefully 
formulated  the  3D  near-plume  problem  using  a  combination  of  heavy  particle  tracking 
and  magnetized  electron  fluid  equations,  and  have  to  a  large  extent  clarified  the  structure 
of  the  solutions,  such  that  numerical  implementation  should  be  greatly  facilitated.  In  area 
(b)  we  have  concentrated  on  an  examination  of  the  far-field  of  the  plume,  where  the 
geomagnetic  field  takes  over  the  dynamics.  Again,  only  a  formulation  and  a  suggested 
method  for  numerical  implementation  were  completed.  Two  graduate  Theses  are  still  in 
progress  in  these  two  areas,  and  we  expect  to  issue  a  supplementary  report  when  they  are 
finalized. 
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1  Near-exit  plume  model 
1.1  Introduction 

A  multiplicity  of  plume  codes  based  on  hybrid  simulation  has  been  developed  by  many 
researchers  [1,2, 3, 4]  these  codes  mature  and  comparisons  are  attempted  to  lab  or  space 
plume  data,  it  becomes  evident  that  one  of  the  essential  ingredients  is  the  distributions  of 
density,  velocity  and  temperature  of  the  ions  at  the  initial  plane,  which  is  either  the 
thruster  exit  plane,  or  a  plane  chosen  some  short  distance  downstream  from  the  exit.  Our 
work  intends  to  improve  present  models  for  initial  plane  distributions  in  Hall  thrusters  by 
solving  accurately  the  near-exit  zone. 


The  near-exit  plume  is  characterized  by  different  aspects  that  are  not  covered  in  the  usual 
plume  models: 

•  Strong  magnetic  field. 

•  Ionization. 

•  3D  effects  in  the  electron  population  induced  by  the  hollow  cathode. 


Such  effects  can  be  taken  into  account  by  a  hybrid  model.  Heavy  particles,  ions  and 
neutrals,  are  modeled  using  Particle-In-Cell  (PIC)  methods,  whereas  electrons  are 
considered  as  a  continuum.  The  main  difference  with  respect  to  previous  hybrid  plume 
model  is  in  the  electron  treatment,  which  includes  the  3D  anisotropy  induced  by  the 
magnetic  field. 


Figure  1.1  Near-exit  plume  diagram.  The  cathode  induces  a  highly  anisotropic  profile.  Due  to  the  fast 
diffusion  along  the  magnetic  field  lines,  its  effect  extends  far. 
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Figure  1.1  presents  the  typical  geometry  in  the  near-exit  plume.  Even  though  the 
geometry  of  the  acceleration  channel  and  the  magnetic  field  is  axi-symmetric,  the  cathode 
breaks  the  symmetry.  The  diffusion  along  the  magnetic  lines  is  much  larger  than  the 
perpendicular  diffusion.  Thus,  the  cathode  imposes  the  electron  temperature  and  the 
electric  potential  on  the  magnetic  lines  that  touch  it.  The  effect  of  the  cathode  extends 
really  far  due  to  the  fast  diffusion  along  the  magnetic  lines.  This  effect  and  its  range  is 
difficult  to  model  because  the  transport  across  the  magnetic  lines  has  two  components 
very  different  in  nature  and  value:  the  diamagnetic  transport,  perpendicular  to  the  driving 
potential  gradient,  modified  by  the  pressure  gradient,  and  the  collisional  perpendicular 
transport.  The  diamagnetic  transport  is  closely  related  to  the  ExB  drift  and  it  is  a 
consequence  of  the  Larmor  motion.  The  collisional  perpendicular  transport,  much  smaller 
in  size,  has  to  do  with  perturbations  on  the  Larmor  motion,  such  as  collisions  or 
turbulence.  Note  that  in  the  axi-symmetric  models,  the  diamagnetic  transport  lays  in  the 
azimuthal  direction,  balancing  out.  However,  in  a  3D  case,  such  as  this  one,  this  transport 
has  to  be  taken  into  account  because  it  is  the  most  important  contribution  in  the  direction 
perpendicular  to  the  magnetic  field. 


In  what  follows  we  discuss  the  electron  fluid  equations  and  the  special  mathematical 
difficulties  that  appear  as  a  consequence  of  the  combination  of  strong  magnetization  and 
three-dimensionality. 


1.2  Electron  equations 

Electrons  are  modeled  as  a  continuum.  The  equations  to  be  solved  are  charge,  momentum 
and  energy  conservation. 

Since  the  Larmor  radius  is  really  small  compared  to  the  typical  length  in  the  problem,  a 
simple  diffusive  model  is  proposed  for  all  these  equations.  Such  an  approximation  is 
acceptable  for  the  motion  perpendicular  to  the  magnetic  field  lines,  but  it  is  not  as  correct 
along  the  lines.  However,  the  proposed  model  is  expected  to  yield  results  that  contain 
most  of  the  physics. 


In  the  next  subsections,  the  different  equations  are  described  and  explained  briefly. 


1.2.1  Charge  conservation 

Assuming  quasineutrality: 

Vj,+Vje=0  (1.1) 

j(.  is  the  ion  current  density,  and  it  is  obtained  in  the  PIC  submodel.  jc  =  -ene\e  is  the 

electron  current  density.  It  is  one  of  the  unknowns  that  must  be  obtained  from  the 
electron  equations. 
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Momentum  conservation 

In  the  diffusive  limit: 

o=-— v(»,r.)+^-v*+i<®j,xb+i»a  (1.2) 

me  me  e  e 

where  <(>  is  the  electrostatic  potential,  Te  the  electron  temperature,  coe  =eB/me  is  the 

gyrofrequency,  6  =  B  /  B  is  the  unit  vector  in  the  same  direction  as  the  magnetic  field  and 
ve  is  a  frequency  that  measures  the  collisionality  and  turbulence. 

can  be  solved  from  this  equation: 

i  =  jj>  +  \eH  +  il  (1-3) 

where 

jeb=-crb b  Vf+-VTe(lnne-l)  (1.4) 

L  e  J 

j,„=-CT„bx  v^+ivr,(lnn,-l)  (1.5) 

L  e  J 

i,i=-^  V^’+iv^Onn.-l)  (1.6) 

L  e  J 

jeb  is  the  current  along  the  magnetic  field  line.  \eH  is  the  diamagnetic  component  of  the 
current.  It  contains  the  electron  current  due  to  the  ExB  drift,  plus  the  curvature  and  VJ? 
drifts.  jei  is  the  collisional  perpendicular  current.  It  accounts  for  collisions  and 
turbulence. 

It  is  important  to  point  out  that  in  axi-symmetric  cases,  jeH  is  aimed  in  the  azimuthal 
direction,  therefore  cancelling  exactly.  Thus,  it  is  not  usually  included  in  2D  models. 

In  these  equations,  the  thermalized potential,  (j>* ,  is  used.  It  is  defined  as: 

0*  =  ^-  — In  ne  (1.7) 

e 

Vx  =  V-66-V  is  the  gradient  in  the  perpendicular  direction  to  the  magnetic  field. 
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The  conductivity  for  the  plasma  is  different  in  the  different  spatial  directions.  In 
equations  1.4, 1.5  and  1.6,  these  conductivities  are  called  ab ,  aH  and  a± : 


e2n„ 


mv 

e  e 


e2ne  coe  _  ene 
me  co ]  +  v2e  B 


(ve  «coe ) 


<7± 


e\  K 

me  co2e  +  v2e 


(1.8) 

(1.9) 

(1-10) 


Considering  that  ve  «  coe  near  the  exit: 

°h  »  o'//  »  oL 

The  ratio  between  these  different  conductivities  is  the  hall  parameter,  fiH  =  coelve.  For 
example,  <Jb  /  aH  ~  ,  or  crH  /  cr±  ~  J3H  .  This  means  that  the  plasma  tends  to 

homogenize  much  faster  along  the  magnetic  lines,  where  ab  is  the  conductivity.  Then, 
along  the  magnetic  lines: 

b- --b- VTe(ln«e  -1)  =  0  (1.11) 


Actually,  from  the  energy  equation,  we  will  obtain  that  6  -  VTe  =  0 ,  which  means  that 
1.11  is: 


6-v^*  =0 


(1.12) 


Even  though  6  -  V<j>'  =0  and  6  •  VTe  =  0 ,  jeb  is  non-zero. 


Plugging  equation  1 .3  into  1.1: 

Vx(6<t„)-v<s-  +v.(<Tiv1f)-vx(6r„).v(r,/e)-v.(rivx(r,/e))  = 

=  v.j,+v-(6/J 

r,  =  <Tj  (In  II,  -1),  where  j  =  H,l. 


It  is  important  to  note  that  in  this  equation  there  are  terms  of  very  different  order  of 
magnitude.  The  terms  that  contain  aH  and  yH  are  larger  by  a  factor  coelve,  provided  the 
other  factors  in  them  are  comparable  to  those  in  other  terms. 


7 


1.2.3  Energy  conservation 

The  energy  conservation  is  given  by: 

d_ 
dt 

where  qe  is  the  electron  heat  conduction  flux,  and  nevjaEi  are  the  ionization  losses. 


(3  ) 

f 

-nT 

+  V- 

V2  e  V 

V 

.  5  T 

-h~— +q, 

2  e 


=  -ie-V(/>-neviaEi 


(1.14) 


The  heat  conduction  can  be  obtained  using  a  diffusive  approximation,  similar  to  the  one 
used  to  obtain  the  electron  current  density  (equation  1.3): 

qe  +  +Qel  (U5) 

where 

<?,s=-v46-vr,  0.16) 

1'H  =-KHf>xVTe  (1.17) 

q  (1-18) 


Thermal  conductivity  is  highly  anisotropic.  The  values  of  the  different  coefficients  are: 


5  neTe 

Kb  =  n 

2meve 

(1.19) 

_  5  neTe  coe 

2  me  co2e  +  v] 

(1.20) 

,  .  _5neTe  ve 

K±  2  me  co2  +  v2e 

(1.21) 

For  o)e  »  ve ,  the  coefficients  have  really  different  orders  of  magnitude: 

Kb  »  Kh  »  K± 

As  it  happened  with  the  electrical  conductivities,  the  ratio  between  the  different  thermal 
conductivities  is  the  Hall  parameter. 


Thermal  diffusion  is  high  in  the  direction  of  the  magnetic  field.  That  means  that,  as  a  first 
approximation,  the  temperature  is  constant  along  the  lines: 

6-V7;=0  (1.22) 
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Even  though  6  •  VTe  is  almost  zero,  qeb  is  not  zero.  That  means  that  equation  1 .14  can  be 
written  as: 


3  d 
—eti  — 
2  e  dt 


+  V- 

bqeh  -  b 

7s  ) 

--In  ne 

T  1 

J  eb 

V  e  ) 

A  2  ) 

e  ) 

+ 


h  1 

T 

- 

—  Inn 

-*—4> 

V2  eJ 

e 

+Vx(bE„)-Vf  +  V-(2±V±f) 

-Vx(br„)-V 


^-(ene)  +  enev\-a^J--^-]nne  -  <f 
dt  V  e  e 


f  t'] 

f 

(T  v 

-V- 

r  v 

1 1 v  i 

k  e  J 

V 

\e  )y 

(1.23) 


The  first  term  in  the  LHS  is  just  the  time  variation  term  for  temperature.  The  second  term 
is  the  energy  transport  along  the  magnetic  field  lines.  We  would  eliminate  this  term  by 
integrating  along  the  magnetic  lines.  Note  that  the  energy  transport  includes  the  potential 
energy,  -ene(f> .  This  potential  energy  transport  can  be  obtained  from  the  Joule  heating 
term,  as  we  will  explain  a  few  paragraphs  below.  The  rest  of  the  LHS  terms  are  transport 
terms  for  energy,  depending  on  (which  induces  electron  current,  and,  thus,  electron 
enthalpy  transport)  and  Te  (which  induces  electron  current  and  heat  conduction).  In  all 

these  terms,  potential  energy  transport  is  also  taken  into  account.  The  RHS  has  two  main 
contributions:  the  energy  variation  due  to  mass  varying  in  time  and  the  energy  change 
associated  to  ionization.  In  the  RHS  the  potential  energy  is  also  accounted  for. 


The  coefficients  in  the  equation  are: 


T  \ 

( 5  ,  1 

2,= 

_ e_\ 

_  e 

V2  7 

r,= 


T 

f  7  1 

-<t>  (ln/ie  —  1) 

—In  ne 

-~ln  ne 

e 

^2  ) 

(1.24) 

(1.25) 


where 


Equation  1 .23  is  obtained  by  using  that: 

=  +  (1.26) 


Taking  into  account  that  charge  conservation  for  the  electron  species  is 

-~(ene)+V-je=-enevn  (1.27) 

dt 

the  final  expression  for  the  Joule  heating  is: 
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(1.28) 


V^  =  -V-(je^)+0— {en^-ferieVi 
ot 


1.2.4  Summary 

Equations  1.12,  1.13,  1.22  and  1.23  can  be  solved  for  <j>' ,  Te ,  jeb  and  qeb. 


1.3  Integration  along  magnetic  lines 

In  this  problem,  two  quantities,  and  Te ,  are  constant  along  the  lines.  Thus,  the 
problem  is,  basically,  2D.  Considering  this,  the  problem  can  be  solved  in  some  other, 
more  natural,  coordinates.  Let  us  define  A, ,  A2  and  5  such  as: 


6- VA,  =0 

(1.29) 

6  •  VA2  =  0 

(1.30) 

—  =  b 

(1.31) 

ds 

where  x  =  (x,  y,  z)  are  the  spatial  coordinates. 


Figure  1.2.  Magnetic  coordinates  A,  and  A2  in  an  axil-symmetric  magnetic  field. 

Note  that  A,  and  A2  are  two  variables  that  define  each  magnetic  field  line  (thus,  these 
variables  are  constant  along  magnetic  field  lines). 


We  require  Vs  -(VA,  x  VA2)  =  0.  We  can  also  impose  VA,  •  VA2  =  0,  but  this  last 
condition  is  very  restrictive  and  can  be  relaxed. 
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As  an  example,  consider  a  magnetic  field  with  cylindrical  symmetry.  In  this  case,  the 
variables  Xt  and  X2  could  be  the  magnetic  flux  function  and  the  azimuthal  angle  (see 
figure  1.2) 


In  this  new  reference  system,  the  divergence  of  a  vector  V  is: 

v.v  =  v,^+v^+vV3V 


ds 


dA, 


dX, 


(1.32) 


In  the  new  reference  system,  the  volume  integrals  must  be  done  considering: 


dx 

dx  dydz  =  dsd\  dX2 - 

ds 


( 


dx  dx 
x 


[dA,  cU2, 


(1.33) 


Then,  the  integral  of  the  divergence  of  V  must  be: 


dxdydzV  •  V  =  dsd) l,  dX: 


dx 

ds 


( dx  av  av  av 

X -  Vs. —  -  +  V/L. —  +  V/L.- 

1  ^  '"i  2 


8  A,  dX2J_  ds 


dX, 


dX 7. 


(1.34) 


To  simplify  this  expression,  first  consider  the  following  equivalences: 


dx 

ds 

dx 

ds 

dx 

ds 


dx  dx 
x 


(SAj  dX2 
f  dx  dx  ^ 


Vs 


Vs 


dx  dx 


f  dx  dx^ 
x 


^  dX2 


"  Vs-(V^xV^) 

~  \VA,xVA2  I 

dX,  dX2 

v^ 

v^ 

dx  ^  dx 

Vs-  (VXi  x  VAj) 

>* 

X 

< 

N> 

K> 

|  ds  dX2 

V  X2 

va2 

dx  dx 

Vs-(VT,xVA2)  \VXixVX2\  ds  dXl 


(1.35) 

(1.36) 

(1.37) 


In  these  equations  we  have  used  VX{  x  VA2  =  6  |  VXt  x  VX2  |  and  b  •  Vs  =  •  Vs  =  1 . 

This  does  not  mean  that  Vs  is  parallel  to  6 . 


Using  equations  1.35,  1.36  and  1.37  to  simplify  1.34: 


dxdydzV  •  V  =  dsdX]  dX2 


d_ 

ds 


f 


Vs- v  d 

+ 


{  Yi 


-V 


JVX]xVX2\J  dX]{\VXlxVX2JJ 

To  get  to  this  equation  we  have  used: 


+  - 


dX, 


vx2-v 

{\VXlxVX2\J 


(1.38) 
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(1.39) 


d 

Vs 

d 

f  VA,  ] 

5 

VA,  ) 

ds 

v.  1 V  A,  x  V  A,  L 

dA, 

ll  VA,  x  VA,  |J 

dA2 

Ll  VA,  x  VA,  |J 

a_ 

ds 


^  dx  dx 
3A,y 


_8_ 

dA, 


dx  dx 
ds  8A2j 


+  ■ 


dA. 


rdx  dx" 
ds  dA, 


2  V 


Equation  1.38  is  important  because  allows  us  to  simplify  the  equations.  Integrating  in  s , 
that  is,  along  magnetic  field  lines: 


f'rfs 

a 

Vs-V  ' 

a 

VA,  •  V  ^ 

a 

'  VA,  •  V  "1 

as 

J  VA,  x  VA2  |v 

_l _ 

aA, 

<1  VA,  x  VA,  L 

1 

aA, 

VA,  x  VA,  |J 

Vs-V  dsB  VA.-V  dsB  VA,-V 


L| VA,  x  VA, 1  dA,  |VA,  x  VA, 1  dA2  | VA,  x  VA2 

a  f  .  ^ 

d 


+ 


+- 


dA,' 


VA,-V 

|  VA,  x  VA,  \) 


dA , 


}  'ds 


va,-v 
I  VA,  x  VA,  \J 


(1-40) 


Here,  sB  (A, ,  A2)  is  the  equation  that  describes  the  boundaries  of  our  domain.  sB  can  be 
s0  or  s, . 


Considering  that  F  =  5  -  sB  (A, ,  A2 )  =  0  is  the  equation  for  the  boundary: 

r)c  a 

VF  =  Vs-— -^-VA,  —  — VA,  =  kN  (1.41) 

dA)  dA2 

where  iSf  is  the  normal  to  the  boundary,  and  k  is  simply  the  magnitude  of  the  vector.  It  is 
also  true  that: 

&b-N  =  b  -  Vs-^-b-  VA,  --^-b-  VA,  =  — •  Vs  =  l  (1.42) 

dAl  '  dA2  2  ds 


Then,  k  =  1  /(6  •  ft)  and  we  can  write: 

ds 


Vs  - 


dA , 


d<>  N 

^VA,-^-VA,=-^ 


dA 7 


b  N 


(1.43) 


Finally,  the  integral  of  the  divergence  of  V  along  the  magnetic  field  lines  can  be  written 
as: 


12 


[  .  dxdydzV  ■  V  =  dA  dA. 

■'alnno  h  ^ 


(b-  N)  |  VA,  x  VA2 


(b  -  N)  |  VA,  x  VA2 


+A  f^_gjL-Y_L_L  J  ,/v  v<  v 

|VA,xVA2|  J  8A2{J‘  o  |VA,xVA2 


(1.44) 


where  Vu  =  V  •  N . 


Let  us  apply  this  to  equation  1.13.  Among  other  terms,  we  will  find: 

jJ  p.i„q„VV(bxVf)l  8  C,±a„VA2-(bxVf) 
aA,(/*  IVA.X.VAJ  J  8A2 Jf»  | VA,  x  VAj  | 


(1.45) 


Considering  that  6  •  V^* 


^  8  A,  n  8  A2 


(1.46) 


We  also  know  that  6  •  (VA,  x  VA2 )  =|  VA,  x  VA2  j .  Then: 

A!  [t‘>fc<TtfV^-(^xV^)l  ,  d(  \s'±a"WX 2-(bxVf)  = 
aA,!/*  |va,xva2|  J  aAAJ*  iva^vaj  , 

a  f—  af  1  a  f—  a^*l  a<rw  a^*  8oh  aa* 

=  -— —  <Tw —  + —  oh-j—  - - —  + - — 

dA\  8A2)  8A2  V  8AI J  8A]  8A2  8A2  8A] 


where 


=  [dsa„ 


(1.47) 


(1.48) 


Another  term  that  will  appear  is: 


|  VAlxVA2\)  8A 


Jo  |VA,xVA2| 


a  f—  af  a  i  —  a^*  a  (—  a^*  a  (—  aa* 

— —  CTx.il— —  +  — —  CTl,12— -  +  —  (71,12 -  H - (7 1,22 - 

O  A  A  n  A  A  A  A  A  /I  A  A  A  A  T  A 


where 


a^v  aAj  aA,^  ’  aAj  8A 


fs,  VA,  -VA, 

i,y  =  J  (Act, - —  for  i  =  1,2,/ =  1,2 

i|VA1xVA2| 


aAj  aA2 


(1.49) 


(1.50) 
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In  a  case  where  VA,  -VA2  =0,the  equations  are  simpler  because  cr±, 12  =  <r±,2i  =0, 
gtj.,11  =  J  350^  |  VA,  |/|  VA2  I  and  01,22  =  j  dscr1  |  VA2 1/|  VA,  |. 


Working  all  terms  similarly,  equation  1.13  becomes: 


Saw  3^*  3o~w  3^*  |  dyH  3  Te  dyH  d  Te  ^ 

3  A,  3A2  3A2  3  A,  3  A,  3A2  e  df  3A,  e 


+- 


_3_ 
'  d\ 


3A, 

3  T 


-  df  -  df  I  3 

(J  1,1 1 - 1-  0  1,12  - — - 1  + 


3A, 


'3AJ  3A2 


df  -  df 

(71,12-7-7“+  <71,22 


\ 


3A, 


-  3  Te 

y±,u  ,  +  /i,i2  -j  * 

3  A,  e  3A2  e 


3A 


Xl.!2 


3A2y 

3 


•2  V 


=-4  f,*  ZAJ 


3 A,  1,  %  |  VA1  x  VAj  |J 

3/W  +  7e/V 


+  - 


_d_f  - _ a 

3 A,  e  ri’22  3A2  e  j 

+  4f"<fc  + 

91;t  '  l^xV^U 

./,W  +  7e5 


+ 


(b-  N)  |  VA,  x  VA2 1 


ii=i0 


(b  •  N)  |  VA,  x  VA2 


(1.51) 


Similarly,  equation  1 .23  becomes: 


p.  3  ene 
2  |  VA,  x  VA2 


35 


44  + 

3/  e 


_j - - - — 


a 


eN 


(b  •  N)  |  VA,  x  VA2  I 


Q, 


eN 


(b  •  N)  |  VA,  x  VA2 


3£tf  3^*  |  3Z/r  df  |  3r//  3  7e  dTH  3  Te 
3A,  3A2  3A2  3 A,  3A,  3A2  e  df  3A,  e 


+- 


3A, 


-  3^*  -  3^ 

2-1,11 - H  2-1,12 


3A, 


3A,  J 


A 


3A, 


-  3^*  -  3^' 

2-1,12  '  +  2-1,22 


♦A 


3A, 


_3_ 

3A, 


f  5  7-  37: 

1 1,11 - 1-1 1,12 


3A,  e 

-I: 


df  e 


3A. 


df  J 
\ 


F  3  7-  37; 

ri,I2- - -+ 11,22— — 

3A,  e  3A2  e  y 


35 


|  VA,  x  VA2 


*-2 


f  ^  \ 

a — --<j) 
\  e  J 


where  QeN  is  the  total  energy  flux  perpendicular  to  the  boundary 


(1.52) 


QeN  ~  QeN  J  eN 


51 

V2  e 


<t> 


(1.53) 
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These  equations  are  equations  for  tp'  and  Te .  Before,  jeb  and  qeb  were  also  unknowns, 
but  integrating  along  the  lines  has  eliminated  these  variables.  This  integration  has  also 
eliminated  the  dependence  with  5  .  Only  the  dependence  with  Xt  and  X2  is  left.  These 
variables  determine  univocally  each  magnetic  line.  The  best  way  to  visualize  these 
variables  is  considering  an  axi-symmetric  magnetic  field.  Then,  any  cylindrical  surface  is 
cut  once  and  only  once  by  each  magnetic  line  and  each  point  in  the  cylindrical  surface 
corresponds  univocally  to  a  magnetic  line.  That  is  what  is  plotted  in  figure  1.2.  Xt  is 
equivalent  to  the  axial  coordinate  in  the  cylinder,  and  X2  is  equivalent  to  the  azimuthal 
angle. 


It  is  interesting  to  take  a  look  into  the  form  of  the  equations.  There  are  terms  derived  from 
the  diamagnetic  transport,  such  that 


8<jh  dtp'  d<7H  dtp' 

d/ij  8X2  8X2  8X, 


(1.54) 


These  terms  are  characterized  by  the  lack  of  second  derivative  and  they  are  usually  the 
higher  order  terms  in  the  equation  (by  a  factor  coe  / ve). 


The  collisional  perpendicular  transport  terms  are  diffusion-like,  such  that 

,  o  1  -  00  -  oc 

+ 


A 

ex, 


d<j>  -  a  A 

CTj.,1 1  — z— -  +  01,12 - 

8X,  dX2 


d  -  dtp'  -  dd) 

-  <71,12  +  (71,22  — 

dX2  v  dX,  dX2  j 


(1.55) 


Since  these  equations  have  been  integrated  along  the  magnetic  lines,  the  flows  across  the 
boundaries  appear,  such  that 

_ JiN  JeN _ 

(b-N)  |  VT,  x  VX2 1 


+  — 


JiN  +  JeN 


(b-N)  |  V2,  xVX2 


(1.56) 


Any  other  term  is  just  the  corresponding  contribution  from  equations  1.13  and  1.23, 
integrated  along  the  magnetic  field  lines. 


1.31  Structure  of  the  solution 


Both  equations  1.51  and  1.52  have  similar  structure.  To  study  the  possible  solutions,  let 
us  take  equation  1.51  in  the  simplified  case  Te  =  const,  and  VT,  .VX2  =  0 : 


dan  d<f)'  8<jh  dip'  8  —  dip' 

- — | - j - <7i,U - 

8Xj  8X2  8X2  8X]  8X,  V  8X]J 


f 


+  ■ 


8X. 


<71,22 


2  V 


dtp' 

8X. 


(1.57) 


2  7 


where  uh,  <7i,u,  <71,22  and  S  are  known  functions  of  X.  and  X2. 
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Looking  at  the  order  of  magnitude  of  the  different  coefficients,  the  terms  that  contain  an 
are  much  larger  than  any  other  term  by  a  factor  coe/ve.  This  means  that,  to  zeroth  order: 

8cjh  dft  dan  dft 

- + - ^-  =  0  (1.58 

dX,  dX2  dX2  8X, 

This  equation  imposes  that  <f> ’  is  constant  along  the  lines  <th  =  constant . 


Figure  1.3.  Hall  thruster  representation  in  a  A,  -  X2  plane.  a  h  —  const,  lines  are  plotted. 

Note  that  this  equation  may  not  (and  usually  will  not)  be  enough  to  determine  (j, T ,  since 
its  variation  across  lines  of  constant  an  is  left  indeterminate  by  the  neglect  of  the  higher 
order  terms.  Let  us  look  at  an  usual  Hall  thruster  configuration  plotted  in  a  A,  -  X2  plane 
(see  figure  1 .3).  Assume  that  A,  and  X2  are  the  axial  and  azimuthal  coordinates  that 
locate  the  magnetic  field  lines  (as  seen  in  figure  1 .2).  A  typical  configuration  of 
oh  =  const,  lines  is  plotted  in  the  A  plane.  The  boundary  conditions  are  imposed  on  the 
upstream  and  downstream  boundaries  A,  =  const,  and  on  the  contour  of  the  image  of  the 
hollow  cathode  in  the  A,  -  X2  plane.  In  addition,  periodicity  is  imposed  on  the  top  and 
bottom  edges.  That  means  that  only  some  of  the  oh  =  const,  lines  have  information  for 
(j) ' .  The  modified  potential  in  the  other  lines  has  to  be  determined  by  the  next  order 
equation,  which  contains  diffusion  terms. 


16 


The  question  about  the  validity  of  this  approximation  naturally  arises.  We  have  to  pay 
special  attention  to  the  boundaries  and  the  points  where  the  gradient  of  gh  is  zero.  In  the 
next  sections  we  describe  the  solutions  at  those  zones. 


Solution  near  the  boundaries 

At  a  non-periodic  boundary  intersected  by  gh  =  const,  lines,  a  thin  layer  may  develop. 
The  physical  meaning  of  such  a  layer  will  be  discussed  later.  Now,  the  mathematical 
solution  in  such  a  layer  will  be  studied. 


gh  =  const. 


Figure  1.4.  Coordinate  system  tangent  to  the  boundaries.  Note  that  the  direction  of  <7 h  =  const. 

lines  is  given  by  angle  a  . 

In  the  2D  space  A,  -  A2  it  is  possible  to  define,  at  least  locally,  two  variables,  £  and  77 
(see  figure  1 .4),  such  that  the  boundary  that  we  are  studying  is  a  line  7  =  const,  and  £  is 
perpendicular  to  7  in  the  sense  that 


^  37" 

- G  1,1 1 

d\  3  A, 


d!;  d?i 
3A2  3A2 


<7i  22  =  0 


(1.59) 


In  such  variables,  and  considering  that  the  gradients  along  7  are  much  larger  than  along 
£ ,  the  equation  for  (j. T  would  be: 


-  if 

01,7717  j  ocy  h 

drj 


(  .  8f  d<j)^ 

sma—  +  cos  or - 

drj  dgj 


=  0 


(1.60) 


where 


37 

Sgh  since  =  - 


—  drj  - 

CF±,1|  +  (71,22 

(1-61) 

V3Aj 

8gh  dr/  i  Sgh  drj 

3A,  3A2  dX2  dXj 

(1.62) 
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(1.63) 


ft 


8a h  cos  a  =  - 


daH  dd, 
dA,  dA2 


+ 


daH  dd, 
dA 2  dA1 


ai,m ,  8a h  and  a  can  be  considered  functions  of  d,  only.  These  parameters  do  not 
change  appreciably  in  t)  -  direction  compared  to  what  they  change  in  d,  -  direction. 


Considering  that  a±,nn  /{d,  8a h)  ~  ve/coe«  1 ,  there  are  too  possible  cases: 

•  a  »  sjve/o)e .  In  this  case,  the  an  =  const,  lines  are  far  from  being  parallel  to  the 
boundary.  That  means  that  the  equation  for  the  solution  near  the  boundary  can  be 
written  as: 


-  d2f  .  df  „ 
ai.rjrj  - - —  +  daH  smfl  — —  =  0 


dr}1  dr} 

Since  01,17,7 ,  8a h  and  a  are  only  functions  of  d, ,  the  general  solution  to  this 
equation  is: 


(1.64) 


a/£\  75/  £\  8a  h  sin  cc 

<p  =  A(g)  +  B(£)ex  p - = - 


7 

V  ) 


(1.65) 


This  solution  only  makes  sense  when  sin  a  >  0 ,  because  otherwise  there  would  be 
an  exponential  growth.  This  means,  basically,  according  to  equation  1 .62: 


dan  dr}  dan  dr}  . 
dA,  dA2  dA2  dA, 


(1.66) 


where  rj(A, ,  A2 )  =  0  is  the  boundary  we  are  interested  in,  and  the  vector 
(dr}/dA,,dr}/dA2)  points  inwards  in  the  plane  A,  -A2. 

If  the  condition  in  equation  1.66  is  not  satisfied,  there  is  continuity  between  the 
boundary  condition  and  the  solution  because  there  is  no  boundary  layer  solution, 
that  is,  the  boundary  layer  solution  has  exponential  growth.  This  means  that  the 
value  of  (j)*  in  each  aH  =  const,  line  is  given  by  the  boundary  conditions  on  the 
boundaries  that  DO  NOT  SATISFY  the  condition  in  equation  1 .66.  The  boundary 
conditions  on  the  boundaries  where  equation  1.66  is  true  do  not  necessarily 
determine  the  thermalized  potential  in  the  crw  =  const,  lines:  a  jump  in  value 
occurs  across  the  thin  layer. 

•  a  ~  -yj ve  l(Oe .  In  this  case,  the  aH  =  const,  lines  are  almost  parallel  to  the 
boundary.  Then,  the  equation  for  <f  is  equation  1.60. 


This  is  a  parabolic  equation,  where  not  only  the  rj  dependence  is  important.  No 
general  conclusions  can  be  drawn  from  this  equation.  However,  the  boundaries 
are  not  expected  in  general  to  be  parallel  to  the  aH  =  const,  lines,  except  in  small 
zones. 
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The  mathematical  solution  of  the  problem  makes  necessary  the  presence  of  these 
boundary  layers.  The  thickness  of  this  layers  is  8~  L  (ve  !coe)  in  general,  and 

8  ~  L^j ve  /coe  when  the  boundary  and  the  <th  =  const,  lines  are  parallel.  L  is  the 
characteristic  length  in  the  problem,  associated  to  the  density  and  magnetic  field 
gradients.  The  thickness  of  the  layers  can  be  written  in  terms  of  the  collision  mean  free 
path,  Amfp ,  and  the  Larmor  radius,  pL .  In  such  variables,  8  ~  pL  {LI Tmfp ) .  If  Amfp is 

assumed  to  be  the  mean  free  path  of  classical  electron-neutral  collisions,  the  thickness  8 
would  be  a  small  fraction  of  the  Larmor  radius,  which  makes  no  physical  sense  since  the 
Larmor  radius  is  the  smallest  characteristic  length  in  the  problem.  To  explain  this 
apparent  contradiction,  we  should  consider  that  the  electron  fluid  equations  as  written  are 
only  valid  if  Amfp  «  L ,  or  if  the  classical  value  of  the  collision  frequency  is  substituted 

by  a  more  accurate  anomalous  transport  coefficient.  In  the  case  of  an  anomalous  transport 
collision  frequency,  the  equivalent  mean  free  path  is  Amfp  ~  20 pL  <  L .  Then,  the 

thickness  of  the  layers  by  the  boundaries  would  be  several  gyroradii,  which  is  an 
acceptable  value.  It  seems  that  these  boundary  layers,  where  the  collisional  perpendicular 
transport  is  the  most  important  contribution,  must  be  related  to  turbulence  induced  with 
wavelength  of  the  order  of  the  Larmor  radius,  since  the  classical  collision  mean  free  path 
cannot  explain  their  presence.  Further  study  is  needed  to  determine  the  true  nature  of 
these  layers. 


Solution  near  extreme  points  and  saddle  points  of  gh 

At  a  point  where  3gh  !3\  =  0  and  3gh  IdXj  =  0,  the  approximation  fails.  The  gradient  of 
gh  can  actually  be  approximated  by: 


3(7  H 

~3A^ 


3oh 

~aIT 


at2 

(1.67) 

at2 

(1.68) 

Locally,  two  new  variables,  £  and  77 ,  can  be  defined  to  simplify  the  equation.  These  new 
variables  are  defined  so  that  the  equation  has  the  following  form: 


—  f 3 2  ft  d2  ft\  rd2 oh3<i ')  d2GHdft 

<7 1  ~ — r-  -) - 7-  —C - - - h  77 - - - —  0 

3rj2J  3%2  3j]  3tj 2  38; 


(1.69) 


In  this  equation,  cr±,  32gh  I38,2  and  32gh  I3t f  are  almost  constant. 


The  solution  of  this  equation  just  shows  a  small  zone  where  the  diffusion  dominates.  This 
zone  is  really  small  and  its  effect  is  negligible. 
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Except  for  the  cases  when  <jh  =  const,  lines  intersect  non-periodic  boundaries,  the 
problem  of  determining  the  variation  across  a h  =  const,  lines  is  still  unsolved.  In 
principle,  integration  along  these  “drift  lines”  will  yield  one-dimensional  differential 
equations  for  the  variation  of  the  mean  variables  in  the  perpendicular  (“diffusion”) 
direction.  This  is  now  being  studied  in  detail,  the  difficulty  being  related  to  the  existence 
of  singular  saddle  points  where  these  directions  are  ill-defined.  This  happens  in  particular 
in  the  vicinity  of  the  cathode’s  image,  because  the  plasma  density,  and  hence  the  quantity 
a H ,  has  there  a  local  maximum,  immersed  in  the  more  general  variations  due  to  the  main 
flow.  Some  other  possible  techniques  are  now  under  consideration,  and  further  results 
will  be  available  when  the  M.S.  thesis  by  F.  Parra  is  completed. 


2  Application  of  fluid  models  to  plumes 

2.1  Introduction 

Work  on  a  fluid  model  aims  to  address  noise  and  statistics  issues  inherent  to  a  discrete 
approach  when  simulating  a  Hall  thruster  plume  expansion.  As  such,  a  simplified  plasma 
model  has  been  chosen  as  an  initial  test  case  to  confirm  that  the  fluid  model  works.  This 
preliminary  case  assumes  quasineutrality,  no  neutrals,  no  magnetic  field,  no  collisions 
and  neglects  electron  inertia.  Applying  these  assumptions,  the  fluid  equations  are  the 
ambipolar  momentum  equation,  plus  the  ion  conservation  equation: 

—  {fnlne\i)+V-(minevi<S>\i+Pe)=0,  (2.1) 

ot 

Pm 

-f  +  V.(„,v,  =0.  (2.2) 

Ot 

After  these  are  solved,  the  potential,  assuming  a  polytropic  model,  follows  from, 

L_ 

Teo 


\ 


(2.3) 


— kTe  -e<f>  =  constant.  (2.4) 

r- 1 


Note  that  Equation  2.3  can  be  regarded  as  an  interpolating  approximation  between 
isentropic  flow  (y  =  5/3  for  electrons)  and  isothermal  flow  {y=  1). 


20 


t 


Because  the  form  of  the  fluid  plasma  equations  is  hyperbolic,  the  discontinuous  Galerkin 
method  has  been  chosen  as  the  preferred  solution  scheme  [5].  The  power  behind 
discontinuous  Galerkin  is  in  the  way  it  merges  advantages  of  both  finite  volume  and 
finite  element  methods.  The  result  affords  higher-order  and  thus  more  accurate  solutions 
at  smaller  computational  cost. 

Solving  the  simplified  plasma  model  axisymmetrically  allows  verification  of  our  fluid 
model  against  existing  analytical  solutions  [6].  In  this  test  case,  the  expansion  of  a 
supersonic  jet  into  vacuum  is  simulated  for  comparison  to  known  asymptotic 
distributions  of  the  plasma  properties.  Figure  2.1  gives  PO  results,  constant  across 
elements,  for  normalized  density  and  Mach  number.  In  these  figures,  a  M~10  jet  is 
issuing  from  the  bottom  left  comer  into  vacuum  conditions.  These  low-order  results  are 
reasonable  —  however,  when  trying  to  attain  higher-order  solutions,  numerical  difficulties 
are  encountered  and  no  solution  is  obtained.  These  problems  remain  to  be  resolved  and 
will  continue  to  be  worked  on. 


x(m)  x(m) 


(a)  Normalized  density 

Figure  2.1.  PO  results. 


(b)  Mach  number 


2.2  Far-field  model 

The  purpose  of  the  far-field  model  is  to  examine  how  the  bulk  plume  evolves  at  greater 
distances  from  the  thruster  exit  where  the  complications  of  the  near-field  have  already 
vanished.  In  this  environment,  different  physics  come  into  play,  most  notably  the  effect 
of  the  geomagnetic  field  on  the  final  configuration  of  the  plume.  The  far-field  model  aims 
to  expand  on  previous  analytical  work  and  create  a  numerical  model  to  study  this  regime. 

Formulation  of  the  equations  for  the  far-field  model  is  similar  to  that  for  the  fluid  model, 
but  does  not  neglect  magnetic  field.  Assumptions  made  are  quasi-neutrality,  static 
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magnetic  field  and  Coulomb-dominated  collision  interactions.  Applying  these,  the 
electron  momentum  equation  becomes, 


VPe  =-ene(E  +  vexB)-meneve\e,  (2.5) 

and  assuming  the  electron  gas  is  perfect,  pe  =  nekTe ,  the  polytropic  relationship  can  be 
used  to  express  VPe  as, 

VPe=-L^nekVTe.  (2.6) 

y- 1 


Using  the  electron  current  density,  je  =  -ene\t  and  E  =  -Vtj> ,  Equation  2.5  can  be 
rearranged  to  yield, 


eny 


</>■ 


7  kTe 
y  - 1  e 


.  „  meve  . 

=  jexB  +  ^je. 


(2.7) 


v  kT 

Defining  </>'  =<f> - -  and  E*  =-V</>\ 

y  - 1  e 


°E*  =  Je  X  ^  +  Je’ 


(2.8) 


where  o  = 


e2n„ 


my. 

e  e 


and  (3  - 


eB 

meve  ' 


In  the  far-field,  only  the  geomagnetic  field  must  be  accounted  for  -  in  general,  the  field 
lines  are  straight  and  their  interaction  with  the  plume  will  be  determined  by  their 
orientation  with  the  thrust  axis.  Assume  for  now  that  BEarth  -  Bz .  Taking  components  of 
Equation  2.8, 


Je 


J  ex 
'  Jey 


j  =  oE * 

J  ez  i 


(2.9) 
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One  can  see  that  anisotropy  in  transport  arises  due  to  the  difference  in  coefficients  of 
terms  parallel  and  perpendicular  to  the  magnetic  field.  Perpendicular  terms  have  an 


additional 


1 


factor.  In  general  (3  »  1 ,  meaning 


1 


1  +  /?2  ’  1 +  J32 

perpendicular  direction  is  much  slower  than  that  in  the  parallel  direction. 


«  1  and  transport  in  the 


Assuming  steady-state,  the  continuity  equation  reduces  to  V  •  je  =  0 .  Plugging  the  current 

density  components  in  and  making  use  of  E*  =  -V^* ,  a  modified  Poisson's  equation  is 
obtained, 
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(2.10) 


The  plan  is  to  approach  the  far-field  model  with  a  hybrid  approach,  treating  ions  as 
particles  and  electrons  as  a  fluid  by  solving  Equation  2.10.  The  ion  particle  distribution  is 
used  to  find  the  plasma  density,  ne .  ne  is  then  input  into  the  electron  equation,  solution 

of  which  provides  the  electric  field  that  allows  integration  of  the  ion  equations  to  provide 
the  ion  distribution  for  the  next  time  step. 


The  primary  difficulty  will  lie  in  solving  Equation  2.10  due  to  the  large  discrepancy  in 
length  scale  terms  caused  by  the  magnetic  field.  The  problem  may  be  made  more 
tractable  by  eliminating  the  anisotropy  through  a  transformation  that  incorporates  the  /? 

term  into  the  z  coordinate,  such  as  —  =  (3  — ,  where  (3  -  f3  +  (3'  {/3  represents  a 

dz'  dz 

constant  (3  over  the  domain  and  (3 '  represents  the  perturbation  from  the  constant  value). 
A  similar  result  would  be  obtained  by  using  a  stretched  grid  with  the  elongation  in  the  z 
coordinate.  Employing  a  coordinate  transformation  or  stretched  grid  alleviates  the  length 
scale  difference  between  x,y  and  z  ,  but  does  not  help  with  terms  involving  f3  and  the 
x  or  y  coordinate.  However,  upon  closer  inspection  of  these  terms  in  Equation  2.10,  it 
can  be  seen  that  they  involve  products  of  Vln(o7/?)  and  V</>’ .  If  <f> *  and  <j! (3  are  related 
to  each  other,  i.e.  (f  =  F(al (3),  these  cross-terms  may  cancel  one  another  out  and  not 
pose  any  difficulties  to  the  numerics.  Such  a  situation  is  plausible  since  both  and  a 
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*  l  t 
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are  inherently  related  to  the  plasma  density,  ne .  The  true  nature  of  the  equation  will  be 
discovered  as  this  portion  of  the  model  is  developed  further. 
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